//	--------------------------------------------------------------------
//	Copyright(C) 2006,2007 Zhao Yukun. All rights reserved.
//	This file is a part of the Xeres system.
//	Xeres 2007
//	Contact : xeres.engine@gmail.com
//	$(reservedInfo)
//	$(reservedDeclaration)
//	--------------------------------------------------------------------
///	@file	<quad.cpp>
///	@path	~/src/lib/math/
///	@date	2008/11/18
///	@brief	.

#include "config/config.h"

#include "lib/math/quad.h"

namespace xeres {

	// Quad
	Quad::Quad( const Vector<float>& p0 , const Vector<float>& p1 , const Vector<float>& p2 , const Vector<float>& p3 )
	{
		m_p[0] = p0;
		m_p[1] = p1;
		m_p[2] = p2;
		m_p[3] = p3;

		// FIXME: verify quad!
	}

	// intersect
	bool Quad::intersect( const Vector<float>& o , const Vector<float>& d , const Vector<float>& ng , bool double_side , _Intersection& in ) const
	{
		return _intersect1( o , d , ng , double_side , in );
	}

	// _intersect1
	bool Quad::_intersect1( const Vector<float>& o , const Vector<float>& d , const Vector<float>& ng , bool double_side , _Intersection& in ) const
	{
		float t , no , nd , qd , a1 , a2 , alpha11 , beta11 , _a1;
		Vector<float> inter;
		float a , u1 , v1 , u2 , v2 , u , v;

		// determine major & minor
		int major , minor;

		int mj = AXIS_X;
		if( fabs( ng[AXIS_Y] ) > fabs( ng[AXIS_X] ) ) mj = AXIS_Y;
		if( fabs( ng[AXIS_Z] ) > fabs( ng[mj] ) ) mj = AXIS_Z;
		switch( mj )
		{
		case AXIS_X:
			major = AXIS_Y;
			minor = AXIS_Z;
			break;
		case AXIS_Y:
			major = AXIS_X;
			minor = AXIS_Z;
			break;
		case AXIS_Z:
			major = AXIS_X;
			minor = AXIS_Y;
			break;
		}

		// N.D
		nd = ng.dot( d );

		// parallel
		if( nd == 0 )
			return false;

		// backside
		if( !double_side && nd > 0 )
			return false;

		// d
		qd = ng.dot( -( m_p[0] ) );

		// n.o
		no = ng.dot( o );

		// t
		t = -( no + qd );
		t = t/nd;

		if( t <= 0 ) return false;
		//if( t >= tmax ) return -1;

		inter = o + d * t;

		// a1
		a1 = area( m_p[0][major] , m_p[0][minor] , m_p[1][major] , m_p[1][minor] , m_p[2][major] , m_p[2][minor] );
		a2 = area( m_p[2][major] , m_p[2][minor] , m_p[3][major] , m_p[3][minor] , m_p[0][major] , m_p[0][minor] );
		alpha11 = area( m_p[1][major] , m_p[1][minor] , m_p[2][major] , m_p[2][minor] , m_p[3][major] , m_p[3][minor] );
		beta11  = area( m_p[0][major] , m_p[0][minor] , m_p[1][major] , m_p[1][minor] , m_p[3][major] , m_p[3][minor] );
		alpha11 /= a1;
		beta11 /= a1;
		_a1 = 1 / a1;

		if( ( a = a1 ) > 0 )
		{
			if( (u1 = area( m_p[1][major] , m_p[1][minor] , m_p[2][major] , m_p[2][minor] , inter[major] , inter[minor] )) < 0 )
				return false;
			if( (v1 = area( m_p[0][major] , m_p[0][minor] , m_p[1][major] , m_p[1][minor] , inter[major] , inter[minor] )) < 0 )
				return false;

			if( (v1 + u1) > a )
				goto NEXT_QUAD_SIDE;
		}
		else
		{
			if( (u1 = area( m_p[1][major] , m_p[1][minor] , m_p[2][major] , m_p[2][minor] , inter[major] , inter[minor] )) > 0 )
				return false;
			if( (v1 = area( m_p[0][major] , m_p[0][minor] , m_p[1][major] , m_p[1][minor] , inter[major] , inter[minor] )) > 0 )
				return false;

			if( (v1 + u1) < a )
				goto NEXT_QUAD_SIDE;
		}
		if( ( a = a2 ) > 0 )
		{
			u2 = area( m_p[2][major] , m_p[2][minor] , m_p[3][major] , m_p[3][minor] , inter[major] , inter[minor] );
			v2 = area( m_p[0][major] , m_p[0][minor] , m_p[2][major] , m_p[2][minor] , inter[major] , inter[minor] );
		}
		else
		{
			u2 = area( m_p[2][major] , m_p[2][minor] , m_p[3][major] , m_p[3][minor] , inter[major] , inter[minor] );
			v2 = area( m_p[0][major] , m_p[0][minor] , m_p[2][major] , m_p[2][minor] , inter[major] , inter[minor] );
		}
		goto QUAD_FOUND;

NEXT_QUAD_SIDE:
		if( ( a = a2 ) > 0 )
		{
			if( (u2 = area( m_p[2][major] , m_p[2][minor] , m_p[3][major] , m_p[3][minor] , inter[major] , inter[minor] )) < 0 )
				return false;
			if( (v2 = area( m_p[0][major] , m_p[0][minor] , m_p[2][major] , m_p[2][minor] , inter[major] , inter[minor] )) < 0 )
				return false;

			if( (v2 + u2) > a )
				return false;
		}
		else
		{
			if( (u2 = area( m_p[2][major] , m_p[2][minor] , m_p[3][major] , m_p[3][minor] , inter[major] , inter[minor] )) > 0 )
				return false;
			if( (v2 = area( m_p[0][major] , m_p[0][minor] , m_p[2][major] , m_p[2][minor] , inter[major] , inter[minor] )) > 0 )
				return false;

			if( (v2 + u2) < a )
				return false;
		}
QUAD_FOUND:

		float alpha = u1 * _a1;
		float beta = v1 * _a1;

		if( is_equal( fabs( alpha11 ) , 1.f ) )
		{
			u = alpha;
			if( is_equal( fabs(beta11) , 1.f ) )
				v = beta;
			else
				v = beta / ( ( u * ( beta11 - 1.f ) ) + 1.f );
		}
		else if( is_equal( fabs(beta11) , 1.f ) )
		{
			v = beta;
			u = alpha / ( ( v * ( alpha11 - 1.f ) ) + 1.f );
		}
		else
		{
			float A = 1.f - beta11;
			float B = ( alpha * ( beta11 - 1.f ) ) - ( beta * ( alpha11 - 1.f ) ) - 1.f;
			float C = alpha;
			float D = ( B * B ) - ( 4.f * A * C );
			float Q = -0.5f * ( B + ( ( B < 0.f ? -1.f : 1.f ) * sqrt( D ) ) );
			u = Q / A;
			if( ( u < 0.f ) || ( u > 1.f ) )
				u = C / Q;
			v = beta / ( ( u * ( beta11 - 1.f ) ) + 1.f );
		}
		clamp_ref<float>( u , 0 , 1 );
		clamp_ref<float>( v , 0 , 1 );

		in.inter = inter;
		in.t = t;
		in.u = u;
		in.v = v;
		return true;
	}

	// _intersect2
	bool Quad::_intersect2( const Vector<float>& o , const Vector<float>& d , const Vector<float>& ng , bool double_side , _Intersection& in ) const
	{
		float t , u , v;
		// Rejects rays that are parallel to Q, and rays that intersect the plane of
		// Q either on the left of the line V00V01 or on the right of the line V00V10.

		Vector<float> E_01 = v10() - v00();
		Vector<float> E_03 = v01() - v00();
		Vector<float> P;
		P.cross( d , E_03 );

		float det = E_01.dot(P);
		if( is_zero( det ) )
			return false;
		float inv_det = float(1.0) / det;
		
		Vector<float> T = o - v00();

		float alpha = T.dot( P ) * inv_det;
		if( alpha < float( 0.0 ) )
			return false;

		// if (alpha > float(1.0)) return false; // Uncomment if VR is used.

		Vector<float> Q;
		Q.cross( T , E_01 );
		float beta = d.dot( Q ) * inv_det;
		if( beta < float(0.0) )
			return false; 

		// if (beta > float(1.0)) return false; // Uncomment if VR is used.

		if( ( alpha + beta ) > float( 1.0 ) )
		{
			// Rejects rays that intersect the plane of Q either on the
			// left of the line V11V10 or on the right of the line V11V01.

			Vector<float> E_23 = v01() - v11();
			Vector<float> E_21 = v10() - v11();
			Vector<float> P_prime;
			P_prime.cross( d , E_21 );
			float det_prime = E_23.dot( P_prime );

			if( is_zero( det_prime ) )
				return false;

			float inv_det_prime = float(1.0) / det_prime;
			Vector<float> T_prime = o - v11();

			float alpha_prime = T_prime.dot( P_prime ) * inv_det_prime;
			
			if( alpha_prime < float(0.0) )
				return false;

			Vector<float> Q_prime;
			Q_prime.cross( T_prime , E_23 );
			float beta_prime = d.dot( Q_prime ) * inv_det_prime;
			
			if( beta_prime < float(0.0) )
				return false;
		}

		// Compute the ray parameter of the intersection point, and
		// reject the ray if it does not hit Q.

		t = E_03.dot( Q ) * inv_det;
		if( t < float(0.0) )
			return false; 

		// Compute the barycentric coordinates of the fourth vertex.
		// These do not depend on the ray, and can be precomputed
		// and stored with the quadrilateral.  

		float alpha_11, beta_11;
		Vector<float> E_02 = v11() - v00();
		Vector<float> n;
		n.cross( E_01, E_03 );

		if( fabs( n.x ) >= fabs( n.y ) && fabs( n.x ) >= fabs( n.z ) )
		{
				alpha_11 = ((E_02.y * E_03.z) - (E_02.z * E_03.y)) / n.x;
				beta_11  = ((E_01.y * E_02.z) - (E_01.z * E_02.y)) / n.x;
		}
		else if( fabs( n.y ) >= fabs( n.x ) && fabs( n.y ) >= fabs( n.z ) )
		{
				alpha_11 = ((E_02.z * E_03.x) - (E_02.x * E_03.z)) / n.y;
				beta_11  = ((E_01.z * E_02.x) - (E_01.x * E_02.z)) / n.y;
		}
		else {

			alpha_11 = ((E_02.x * E_03.y) - (E_02.y * E_03.x)) / n.z;
			beta_11  = ((E_01.x * E_02.y) - (E_01.y * E_02.x)) / n.z;
		}

		// Compute the bilinear coordinates of the intersection point.

		if( is_equal( alpha_11 , float(1.0) ) )
		{
			// Q is a trapezium.
			u = alpha;
			if( is_equal( beta_11 , float(1.0) ) )
				v = beta; // Q is a parallelogram.
			else
				v = beta / ((u * (beta_11 - float(1.0))) + float(1.0)); // Q is a trapezium.
		}
		else if( is_equal( beta_11 , float(1.0) ) )
		{
			// Q is a trapezium.
			v = beta;
			u = alpha / ((v * (alpha_11 - float(1.0))) + float(1.0));
		}
		else {

			float A = float(1.0) - beta_11;
			float B = (alpha * (beta_11 - float(1.0)))
				- (beta * (alpha_11 - float(1.0))) - float(1.0);
			float C = alpha;
			float D = (B * B) - (float(4.0) * A * C);
			float Q = float(-0.5) * (B + ((B < float(0.0) ? float(-1.0) : float(1.0)) * sqrtf(D)));
			u = Q / A;
			if ((u < float(0.0)) || (u > float(1.0))) u = C / Q;
			v = beta / ((u * (beta_11 - float(1.0))) + float(1.0)); 
		}

		in.inter = o + d * t;
		in.t = t;
		in.u = u;
		in.v = v;
		return true;
	}

} // namespace xeres
